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Preface 


USAFETAC  prepared  thin  report  to  document  the  Solar  Illumination  Calculator 
(SIC)  computer  program  developed  by  USAFETAC  and  eurrontly  being  used  by  USAFETAC 
and  the  Prototype  World  Wide  Military  Command  and  Control  System  (WWMCCS)  Inter- 
computer  Network  (PWIN).  Thla  report  documents  the  equations  used  by  the  SIC. 

In  the  event  that  this  report  is  incorporated  into  another  report  by  any  agency, 
request  that  USAPETAC  be  furnished  a copy  of  the  new  report  in  all  cases  where  such 
dissemination  is  nbt  prohibited. 

Department  of  Defense  agencies  and/or  their  contractors  who  require  more 
information  about  this  report  should  contact  USAPETAC  for  consultation. 


SOLAR  ILLUMINATION  CALCULATOR 


Introduction 

111 i 

The  basic  assumption  used  in  writing  the  Solar  Illumination  Calculator  (SIC) 
was  that  sunrise,  sunset,  and  twilight  times  could  be  accurately  calculated  for  any 
point  on  the  surface  of  the  earth  if  the  latitude  and  longitude  of  that  point  and 
the  solar  declination  for  the  given  day  were  known.  Work  commenced  on  the  SIC  in 
late  February  1976.  It  soon  became  evident  that  the  SIC  was  a much  more  complex 
project  than  was  originally  envisioned.  This  explanation  will  aid  the  user  in  under 
standing  the  limitations  of  the  SIC,  how  to  use  it,  and  answer  some  questions  as  to 
how  the  equations  were  derived. 

Solar  Declination 

In  order  to  make  the  SIC  more  flexible,  the  author  decided  to  model  the  sun's 
apparent  daily  change  in  declination,  thus  removing  the  problem  of  creating  and 
maintaining  a large  data  file.  It  was  thought,  at  first,  that  the  declination  of 
the  sun,  with  respect  to  a viewer  located  on  the  earth's  surface,  could  be  approxi- 
mated by  a simple  harmonic  equation  of  the  form 


y - A Sin  (2tt  t)  (1) 


where  A represents  the  maximum  declination  of  the  sun  and  T is  a time  factor  express 
ing  the  fractional  value  of  orbital  completion.  The  value  of  A,  called  the  obli- 
quity of  the  ecliptic,  can  be  obtained  from  the  equation 


A * 23027'8'.'26  - 0V4684  (t-1900)  (2) 


where  t=>  the  year.  A natural  starting  point  for  T is  the  March  equinox.  There- 
fore, T can  be  expressed  as  the  time  (in  days)  since  the  last  March  equinox  divided 
by  the  length  of  one  tropical  year  (3o5-24  days).  The  apparent  daily  solar  decli- 
nation can  be  expressed  as: 


DECLINATION  = A Sin 


To  ^ 
365- 24 J 


(3) 


There  are  always  365*24  days  between  March  equinoxes.  When  a leap  year  occurs,  the 
March  equinox  occurs  one  (actually  closer  to  0.75)  day  earlier  than  on  the  previous 
year.  For  example,  the  March  equinox  for  1975  occurred  at  0600Z  on  21  March,  while 
the  March  equinox  for  1976  occurred  at  1150Z  on  20  March.  Note  that  the  time  of 
occurrence  of  the  next  March  equinox  can  be  approximated  by  adding  0.2422  day  to 
the  last  March  equinox  and  subtracting  one  if  the  year  is  a leap  year. 

Compared  with  data  extracted  from  The  Nautical  Almanac  (6),  the  declination 
values  obtained  with  Equation  (3)  are  close,  but  not  exact.  The  majeu^cause  of 
error  is  the  varying  orbital  velocity  of  the  earth.  Kepler's  laws  of  planetary 
motions  explain  this  effect.  However,  if  the  earth's  orbit  is  viewed  as  consisting 
of  four  seasons  of  varying  lengths,  with  the  length  of  each  season  fixed,  a set  of 
four  simple  harmonic  equations  can  be  derived. 


DECLINATION  • A Sir.  (J 


for  period  from  March  equinox  to  June  solstice: 


T 

DECLINATION  = A Cos  (j  ^MEr) 


(5) 


for  period  from  June  solstice  to  September  equinox; 

T 

DECLINATION  = - A Sin  (j 

for  period  from  September  equinox  to  December  solstice; 


DECLINATION  = A Sin 


To  \ 


(6) 


(7) 


for  period  from  December  solstice  to  March  equinox; 

where  T^  = time  since  March  equinox 

SPRING  = 92.78  (number  of  days  in  Northern  Hemisphere  spring  season) 

Ts  = TQ  - SPRING  (number  of  days  since  summer  started) 

SUMMER  = 93*64  (number  of  days  in  Northern  Hemisphere  summer  season) 

Tf  = To  - (summer  and  spring)  (number  of  days  since  autumn  started) 

AUTUMN  = 89.83  (number  of  days  in  Northern  Hemisphere  autumn) 

The  maximum  error  found  in  checking  the  results  of  these  equations  against  data  for 
the  1975-1976  time  frame  is  32  minutes  of  arc.  Expected  error  rates  are  14.8, 

3.0,  17.7,  and  7.7  minutes  of  arc  for  the  spring,  summer,  autumn,  and  winter  sea- 
sons, respectively. 

Equation  of  Time 

An  additional  factor  that  must  be  considered  is  the  equation  of  time  (EOT). 
Basically  stated,  the  EOT  is  the  difference  between  apparent  time  and  mean  time. 

An  in-depth  explanation  of  this  phenomenon  can  be  found  in  any  introduction  to 
astronomy  text  such  as  that  by  Russell  [3J.  The  EOT  results  from  two  effects;  the 
obliquity  of  the  ecliptic,  and  the  eccentricity  of  the  earth's  orbit. 


Klgur*-  I.  Tb*'  Equation  ■'!'  Time 
( !'  -.-'nn  I.  L | l j p . l4y  ) , 


The  maximum  effect  of  the  obliquity  term  is  10  minutes  of  time,  and  of  the  eccen- 
tricity term,  7-3/4  minutes  of  time.  Note  that  both  terms  approach  zero  at  about 
the  timee  of  the  solstices  (for  the  eccentricity  term,  it  is  actually  perigee  and 
aphelion,  which  is  on  the  order  of  11  days  after  the  solstice).  Using  the  winter 
solstice  and  perigee  as  a starting  point  in  time,  two  simple  harmonic  equations  can 
be  derived,  which,  when  combined,  mathematically  measure  this  phenomenon.  Thtv 
NOT  reduces  to  the  forms 


EOT  ■ 10  Sin 


- 11) 

■ tff"111 

Twb 


1 


where  T ■ time  (in  days  since  laet  winter  solstice,  and  182.62  ■ $ tropical  year. 

0 Wfl 

Results  from  this  equation  are  in  close  agreement  with  data  extracted  from  The 
Nautical  Almanac  (6).  A maximum  error  of  1 minute  of  time  occurs  in  the  March  and 
the  late' ' September  through  early  October  time  periods.  Soule  [4]  argued  that  the 
EOT  at  both  sunrise  and  sunset  should  be  computed  to  improve  accuracy.  USAFETAC 
adopted  this  convention. 

Semiduration 

The  semiduration  of  a phenomenon  in  terms  of  the  cosine  of  the  Oreenwioh  Hour 
Angle  (oHA)  is  computed  using  a simplified  version  of  the  time-sight  formula  as 
presented  by  Bowditch  [1]. 


Cos  (OHA)  ■ | - cos  (tatltudefn*^os  (declination)] 

- Tan  (Latitude)  * Tan  (Declination)  (9) 


Alpha  has  a valuo  of  0.833  degrees  for  sunriao/sunset,  6 dogreas,  Id  degroen,  and 
18  decrees  for  civil,  nautical,  and  astronomical  twilight,  respectively,  when  com- 
puting data  at  sea  level  and  with  no  terrain  effect,  These  are  standard,  defined 
values.  Tho  value,  0.833  degrees,  Is  arrived  at  by  adding  the  solar  diameter, 
assumed  to  l.e  fixed  at  16  minutes,  to  an  altitude  correction  factor  for  refraction 
of  34  minutes.  No  methods  are  available  to  account  for  a changing  atmospheric  re- 
fractive index  caused  by  temperature  and/or  pressure  changes.  The  value  of  alpha 
can  be  changed  for  Increases  in  altitude.  This  correction  is  made  solely  for 
sunriBn/sunsot  and  civil  twilight  calculations.  The  equations  UBOd  to  compute  the 
new  valu-’  of  alpha  are  based  on  data  extracted  from  The  Air  Almanac  [5]. 

In  computing  semiduration  there  are  times  when  the  Cosine  (OHA)  exceeds  * 1. 

When  this  happens,  it  indicates  nonoocurrence  of  the  phenomenon,  f nBider  the 
following  examples: 

• Cos  (OHA)  " 1.1  phenomenon  ■ sunriDe 

This  indicates  that  the  sun  will  not  riBO,  because  the  sun  is  always  below 
the  horizon  (winter  cane), 

• Con  (OHA)  ■ - 1,1  phenomenon  * ounrlBe 

This  indicate b that  the  sun  will  not  rise,  because  the  sun  1b  always  above 
the  horizon.  In  addition,  the  cosine  values  for  civil,  nautical,  and  astro- 
nomical twilight  will  also  be  less  than  -1.1  (absolute  value  will  be  greater 
than  J.l).  But,  since  the  sun  doesn't  set,  these  phenomena  will  also  not 
occur  (rummer  case). 

In  the  more  general  Bcnse,  there  arc  two  possibilities  when  computing  oemidurnti on; 
tho  possibility  that  the  phenomenon  doeB  not  occur,  ami  the  poBtilb.il.lty  that  tin* 
phenomenon  occurs  at  leant  once  during  the-  day.  In  the  former,  there  are  two  can's, 
the  winter  case  and  the  nummer  case.  In  tho  latter,  there  are  three  canes j the 


normal  case,  the  spring  case,  and  the  autumn  case.  All  oases  are  based  on  Northern 
Hemisphere  seasons.  In  the  winter  case,  the  Cos  (QHA)  value  is  always  greater  than 
1.0.  The  sun  is  always  below  the  horizon,  hence  no  sunrise  or  sunset,  and  the* 
twilights  may  or  may  not  occur.  In  the  summer  case,  the  Cos  (CKA)  value  is  always 
Iobb  than  "1,0.  The  sun  is  always  above  the  horizon,  hence  no  eunriBe  or  sunset, 
and  the  twilightB  never  occur.  In  the  normal  case,  the  Cos  (QHA)  value  always  lies 
in  the  range  of  equal  to  or  greater  than  -1.0  and  equal  to  or  less  than  1.0,  hence 
the  sun  rises  and  sets,  and  the  twilights  occur.  In  the  spring  case  the  Cos  (UHA) 
value  for  the  beginning  time  of  some  phenomenon  Is  always  In  the  range  of  values 
for  the  normal  caBe,  but  the  Cos  (QHA)  value  for  the  ending  time  is  Ibbb  than  -1.0 
(the  phenomenon  begins  but  does  not  end  before  the  day  is  over).  The  autumn  case 
Is  the  reverse  of  the  spring  case  (the  phenomenon  is  occurring  at  OOOOZ,  but  enda 
before  2h00z)* 

Half  Phenomenon  Length 

Once  the  semiduration  is  computed,  the  SIC  program  computes  half  phenomenon 
length  (HDL) . If  the  Coe  (QHA)  has  an  absolute  value  greater  than  1.0.  the  HDL  is 
sat  to  zero  and  a flag  is  Bet  in  the  program.  If  the  Cob  (QHA)  value  Hbb  in  the 
normal  case  range,  the  aroos  of  the  QHA  is  taken.  ThiB  value  ie  then  multiplied 
by  57.2958  to  convert  a radian  value  to  degrees  and  then  divided  by  15  to  convert 
a degree  value  to  an  xx.yyyy  hour  format.  At  this  point  HDL  has  a value  greater 
than  0.0,  but  less  than  12,0. 

Starting  and  Ending  Timea 

Phenomena  starting  and  ending  times  can  now  be  computed.  If  the  HDI.  has  a 
value  of  0.,  no  computation  is  done.  If  the  HDL  has  a non-zero  value,  starting 
and/or  ending  phenomena  times  are  computed  using  the  following  formulas i 


ST  - 12.  + EOTR  - HDL  (I)  + ADJUST 


(10) 


where  ST  ■ starting  time 
12.  “ local  noon 

EOTR  ■ equation  of  time  at  sunrise 

HDL  (I)  “ half  phenomenon  length  starting 

ADJUST  ■ time  adjustment  factor  for  being  East  or  West  of  areenwtoh  Meridian 
(*  Longitude/15) 

ET  - 12.  + EOTS  + HDL  (I)  + ADJUST  (11) 

whore  ET  ■ ending  time 
12.  » local  noun 

HOTS  ■ equation  of  time  at  sunset 
HDL  (1)  half  phenomenon  length  ending  time 
ADJUST  ■ name  as  above 

Table  1 cuntulrin  the  maximum  error  in  minutes  of  time  for  values  computed  at  ilf'cn- 
wleh.  Since  all  times  arc  based  on  OMP,  ST  and  KT  values  lie  in  the  rang*’  ol’  -Vh , 0 
to  +dB.O.  A value  of  less  than  0.  indicates  that  the  phenomenon  startc/endo  on 
the  previous  Greenwich  Day.  while  a value  of  greater  than  2H.0  indicates  that  the 
phenomenon  utarta/endo  on  the  following  Qreenwlch  Day.  Table  2 contains  an  example. 


Table  1.  Maximum  Timinft  Errors  at  Greenwich.  (Unless  otherwise 


stated, 

equator, 

all  latitudes  are  given  in  degree n 
SUNRISE/SUNSET  AND  CIVIL  TWILIGHT 

N and  s of  the 

NAUTICAL  TWILIGHT 

JAN 

1 

minute  to  40® 

2 

minutes  above  4o° 

3 

minutes 

above 

40' 

m 

1 

minute  to  40° 

2 

minutes  above  40 

4 

minutes 

above 

40' 

MAR 

1 

minute  to  50° 

3 

minutes  above  50° 

4 

minutes 

above 

50' 

APR 

1 

minute  to  30° 

2 

minutes  to  40® 

4 

minutes  to  50* 

6 

minutes  to  70® 

8 

minutes 

above 

50' 

MAY 

1 

minute  to  4o® 

2 

minutes  above  40® 

17 

minutes  at  70* 

8 

minutes 

above 

50' 

JUN 

1 

minute  to  40® 

2 

minutes  above  40® 

2 

minutes 

JUL 

1 

minute 

1 

minute 

AUG 

1 

minute  to  60* 

3 

minutes  at  70® 

8 

minute b 

at  70' 

b 

SEP 

1 

minute  to  60* 

2 

minutes  below  30®s 

3 

minutes  at  60*s  and 

at  ?0*N 

22 

minutes 

at  ?0 

P 

OCT 

1 

minute  to  20® 

2 

minutes  above  20“ 

5 

minute s 

above 

20' 

NOV 

1 

minute  to  30® 

2 

minutes  to  40° 

minutes  above  40* 

5 

minute n 

above 

40 

DEC 

1 

minute  to  50® 

n 

minutes  above  50® 

4 

minuter 

above 

50 

NOTH:  Difference  between  computed  Sic  valuoe  ancl  data  extracted  from  Thu  Nautical 

Almanac  (fj  urine  1975  and  197<>  as  base  years.  Data  were  chncKed  I'rom  H3*T!  to 
70*N  in  10  Inc  remanta.  Notice  the  valueB  for  April  and  November  when  the 
equation  of  time  wog  off  by  1 minute.  The  large  error  in  nautical  twilight 
during  Suptamber  was  due  to  an  autumn  case  occurrence  — the  CTC  produced  a 
normal  cane,  i.u.,  starting  time  of  0022  GMT. 


Table  2.  Example  SIC  Results  for  5 January  1975. 

TIME  (QMT) 

20°N  LAT/120*g  LONO  P0°N  LAT/120*W  LONG 


Astronomical  Twilight 

-0242  (42110) 

1310  (51318) 

Nautical  Twilight 

-0215  (42145) 

1348  (61345) 

Civil  Twilight. 

-0148  (42212) 

1412  (51412) 

Sunrire 

-0124  (42336) 

1436  (51436) 

Sunset 

0934  ('  0934) 

2534  (60134) 

Civil  Twilight 

0950  (50958) 

2658  (60158 ) 

Nautical  Twilight 

10P()  (51026) 

26 26  (60226 ) 

Astronomical  Twilight 

1093  (51053) 

2653  (60253) 

NOTE:  The  numbers  in  parentheses  are  In  a DDHHMM  (day,  hour,  minute)  format,  and  they 

are  printed  out  by  the  printing  routine  of  tha  SIC.  Therefore,  at  20°N  lati- 
tude, 120®W  longitude  B’nset  occurred  at  0134  OWT  on  f>  January  (1734  local  on 
5 January), 


Elevation  and  Azimuth 


Elevation  and  azimuth  angles  are  computed  using  the  following  equations  (from 
Soule  [4]): 


ELEVAT  = ARSIN  [SIN  (STALAT)  x SIN  (DEC) 

+ COS  (STALAT)  x COS  (DEC)  X DOS  (LHA))  (12) 


where  ELEVAT  = Elevation 

STALAT  = Station  latitude 
DEC  = Declination 

LHA  = Local  Hour  Angle  (assumed  to  be  LHA  at  sunrise,  unless  specified 
otherwise  through  parameter  HOUR) 

TSIN  (STALAT)  x SIN  (LHA) 1 rich 

AZIMUT  = ARSIN  [_ * a6g  (fiLEVATT  J K 

where  AZIMUT  = Azimuth  and  STA..AT.  LHA,  and  ELEVAT  are  the  same  as  above.  (If  LHA 
at  sunrise  is  used,  ELEVAT  = 0.0.)  Elevation  and  azimuth  angles  are  not  computed 
for  the  winter  case  (when  the  sun  is  never  above  the  horizon),  or  for  hours  that 
are  either  earlier  than  sunrise  or  later  than  sunset.  Computer-derived  values  for 
azimuth  are  expressed  in  terms  of  the  principle  sine  value.  In  order  to  maintain 
standard  references:  North -i-s-  -zero  degrees;  East  is  90  degrees;  South  is  180 

degrees;  West  is  270  degrees.  The  following  equations  were  taken  from  Dave,  Hal- 
pern,  and  Ityers  [2]. 


AZIMUTH  » PI  - AZIMUT  (14) 


if  the  condition  Cos  (LHA)  < Tan  (DEC) /Tan  (STALAT)  is  true.  Otherwise 


AZIMUTH  = 2xPI  - AZIMUT  (15) 

where  LHA,  DEC,  STALAT,  and  AZIMUT  are  the  same  as  above,  PI  - 3-l4l6,  and  AZIMUTH 
is  true  AZIMUTH.  (In  the  subroutines  no  distinction  is  made  between  AZIMUTH  and 
AZIMUT.)  When  the  SIC  ran  on  the  MAC  Information  Management  System  Honeywell  6o8o 
computer,  library  differences  between  IBM  and  Honeywell  resulted  in  the  sun  rising 
in  the  west  and  setting  in  the  east.  Therefore,  the  azimuth  Equations  (14)  and  (15) 
were  changed  to: 

AZIMUTH  = PI  - AZIMUT  (16) 


if  STALAT/'DEC  is  greater  than  1.  Otherwise 


AZIMUTH  = AZIMUT  (17) 


The  sunrise  azimuth  angle  results  are  in  close  agreement  with  values  computed 
using  a Hewlett-Packard  65  calculator.  (The  H-P  65  azimuth  values  agree  to  within 
± 0.01  degrees  with  those  calculated  by  the  National  Severe  Storms  Laboratory 
celestial  positions  program.)  At  present  the  USAFETAC  has  no  method  for  verifying 
the  accuracy  of  the  elevation  angle  or  the  azimuth  angle  at  some  time  other  than 
sunrise/sunset.  (While  sunset  is  not  computed  directly  it  can  be  found  by  sub- 


tracting  the  sunrise  azimuth  value  from  360  degrees.)  It  should  also  be  noted  that 
solar  declination  is  assumed  to  be  a constant  in  computing  the  elevation  angle. 

This  assumption  is  valid  because  the  maximum  rate  of  change  for  solar  declination 
is  1 minute  of  arc  per  hour. 
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